#!perl
use File::Basename;
die "perl blast100k.pl <fasta> <outprefix> <thread> <blastall> <infofile>\n" if(@ARGV != 5);

open FA, $ARGV[0];
$name = basename($ARGV[1]);
$dir = dirname($ARGV[1]);
$info = $ARGV[4];
$flag = 0;
while(<FA>)
{
	if(/>/){
		if($flag == 0){
			$n+=1;
		}
		if($n==100001 && $ARGV[3] eq "false"){
			$flag = 1;
			last;
		}
		if($flag == 0){
			push @fasta, $_;
		}
	}else{
		if($flag == 0){
			push @fasta, $_;
		}
	}
}
$num = int($n / 1000);
open FB, ">$ARGV[1].${num}k.blastNT.fasta";
foreach $i(@fasta){
	print FB $i;
}
open FC, ">$ARGV[1].${num}k.blastNT.sh";
print FC "cd $dir
blastn -task megablast -query $name.${num}k.blastNT.fasta -db /data/reference/pathoseq/ncbi_nt_v5/20201007/fa/nt.remove.chordata.embryophyta.classification.genome -num_alignments 10 -max_hsps 1 -num_threads $ARGV[2] -evalue 1e-10 -outfmt '6 qseqid sseqid pident length mismatch gapopen qlen qstart qend slen sstart send evalue bitscore staxid ssciname scomname stitle sblastname sskingdom' -out $name.${num}k.nt_v5.blast.txt
#cut -f 1,2,3,7,8,9,11,12,18 $name.${num}k.nt_v5.blast.txt > $name.${num}k.checkblast.txt
#perl /data/softwares/mngs_scripts/blast2nwk.pl --inBlastFile $name.${num}k.nt_v5.blast.txt > $name.${num}k.nt_v5.parse.xls
/home/yong_sun/bin/workflow/venv/bin/python /home/yong_sun/bin/workflow/parseBlast.py $name.${num}k.nt_v5.blast.txt
/home/yong_sun/bin/workflow/venv/bin/python /home/yong_sun/bin/workflow/cal_coverage.py $name.${num}k.nt_v5.uniqblast
/home/yong_sun/bin/workflow/venv/bin/python /home/yong_sun/bin/workflow/calcNTRelAbund.py -nt $name.${num}k.nt_v5.parse.xls -info $info --cofile ${name}_coverage.xls -o $name.nt.ra.xls \n";
